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Abstract 

In  this  work,  the  station  keeping  parameters  at  the 
earth-moon  libration  point,  L4,  were  studied.  First,  the 
equations  of  motion  for  a three-dimensional,  four  body  system 
with  elliptical  orbits  were  derived.  These  equations  were 
then  linearized  about  the  lA  point;  and  optimaLL  control  theory 
was  applied  to  obtain  a linear  feedback  controller. 

The  major  computations  of  the  controller  were  associated 
with  the  gain  matrix,  which  is  the  solution  to  the  time  vary- 
ing Riccati  equation.  Because  of  the  periodic  nature  of  the 
time  varying  gains,  it  was  felt  that  a modified  (fixed  gain) 
control  could  be  used.  The  modified  controller  was  found 
by  computing  the  steady-state  average  of  the  time  varying 
gains. 

Several  observations  were  made  in  studying  the  perform- 
ance of  the  satellite  in  the  vicinity  of  the  L4  point.  First, 
it  was  found  that  the  modified  controller  was  computationally 
much  simpler  than  the  optimum  controller  while  providing  near 
optimail  performance.  Second,  there  is  approximately  a linear 
relationship,  up  to  a point,  between  stationkeeping  cost  and 
distance  from  the  L4  point.  Third,  there  are  initial  solar 
system  configurations  which  minimize  station  keeping  costs. 
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STATION  KEEPING  AT  THE  L4  LIBRATION 
POINT I A THREE-DIMENSIONAL  STUDY 


I.  Introduction 

The  lunar  libration  points  are  equilibrium  positions 
where  a particle,  if  given  the  proper  initial  velocity,  will 
remain  fixed  relative  to  the  other  two  masses  in  a restricted 
three  body  system  (see  Figure  1).  The  collinear  libration 
points,  L^,  Lg,  Lj,  lie  along  the  earth-moon  line  and  have 
been  shown  to  be  vinstable  equilibriiam  points.  The  equilateral 
triangle  points,  Li^,  L^,  have  sides  equeil  to  the  instanta- 
neous earth-moon  distance  and  are  stable  equilibrium  points 
for  mass  ratios  of  less  than  approximately  .04.  While  Lj^  is 
the  primary  point  studied  in  this  work,  the  results  could  be  1 

easily  applied  to  by  symmetry.  | 

In -the  solar  system  there  are  forces  acting  upon  the 
satellite  other  than  just  the  earth  and  moon.  The  libration 
points  are  no  longer  equilibrium  positions.  It  is,  however, 
convenient  to  think  about  the  positions  in  space  that  corres- 
pond to  the  libration  points.  Stability  and  station  keeping 
can  be  studied  with  reference  to  this  point. 

The  lunar  libration  points,  and  in  particular,  have  i 

potential  value  as  the  exploration  of  space  is  expanded.  In 
an  Air  War  College  report,  Lt  Col  Berge  (Ref  2:45-65)  con- 
firmed that  the  point  has  value  as  a communications  link 
for  interplanetary  missions,  back-side  lunar  missions,  and  a 

l! 
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Figure  1 

Restricted  Three  Body  Libration  Points 
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variety  of  scientific  missions. 

The  object  of  this  thesis  is  to  continue  the  work  of 
previous  AFIT  theses  (Ref  6j  7j  8;  10) . The  equations  of 
motion  with  respect  to  the  Ii4  libration  point  will  be  derived 
based  upon  a three-dimensional  four  body  system  (earth-moon- 
svin-satellite)  with  elliptical  orbits.  And  station  keeping 
parameters  will  be  examined  relative  to  this  model. 


FormxAlation  of  the  Equations  of  Motion 

Considering  the  system  of  four  bodies  of  Figure  2,  the 
equations  of  motion  referenced  to  the  inertial,  sun-centered 
coordinate  system  (I,J,K)  can  be  written  as 


(1) 


where  m_  is  the  mass  of  the  satellite 
P 

r _ is  the  vector  from  the  sun  to  the  satellite 
sp 

Fg,  Fjyj,  Fg  are  the  gravitational  forces  of  the  earth, 
moon  and  sun  referenced  to  the  inertial  frame 
The  equations  of  motion  can  be  written  in  terms  of  a 
rotating  coordinate  system  (x,y,z)  centered  at  the  earth-moon 
barycenter  (Ref  Figure  2)  by  using  the  transformation  matrix 
(Ref  13*432-435) 


where 


CoS  0 CoS(9  + 0)  CoS  & Sin(©  + 0)  Sin  3 
-Sin  (9  + 0)  CoS  (0  + j2()  0 

-Sin  3 CoS(e  + 0)  -Sin  3 Sin(e  + 0)  CoS  3 


3 is  the  inclination  of  the  earth-moon  plane  to  the 
ecliptic 
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e is  the  true  anomaly  of  the  earth-lunar  orbit  about 
the  sun 

0 is  the  svin  direction  relative  to  rotating  frame 
In  terms  of  the  rotating  coordinate  system,  the  equations  of 
motion  are 


A r„^  = + P + F 

p - sp  e m s 


where  the  gravitational  forces  referenced  to  the  rotating 


coordinate  system  are 

p = - Is  K , F = _ J 

“^e  2 ^ep*  'm  3 

ep  ^mp 


r « F = - r 

^mp  ' ^s  -A--  P ^sp 

r-^ 

sp 


and  (Ref  Figure  3) 


where 


= (Tq  + €)  + + p^j  r^  = a^  (1  - e Cos  ji) 


r is  the  earth  - barycenter  distance 

r^jj  is  the  moon  - beirycenter  distance 

r^^  is  the  earth  - satellite  distance 
ep 

r„„  is  the  moon  - satellite  distance 
mp 

a^  is  the  semi-major  axis  of  the  earth  orbit 
Sjjj  is  the  semi-major  axis  of  the  lunar  orbit 
e lunar  eccentricity 
U - lunar  eccentric  anomady 

6,t|f,o-  are  the  components  of  the  vector  from  earth- 
moon  barycenter  to  satellite,  R,  in  the  x,y,z-frame 
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Butt  from  Figure  2 


r„_  = r,  + r 
sp  80 


Differentiating 


*^8b  “ - **80 


r = (u)^  + ^ 


where  the  cross-product  matrices  (i)_  and  S)  are  defined  as 
(Ref  13«20) 


0 

• 

-e 

0 

0 

e e 

-0 

0 

« 

e 

0 

0 

• t • 

I (0  = 0 

* X 

0 

0 

0 

0 

0 

0 

0 

0 

where 


r„^  - vector  from  sun  to  satellite 
sp 


- vector  from  svin  to  earth-moon  barycenter 

referenced  to  rotating  coordinate  system  (x,y,z) 
Tg^  - vector  from  svm  to  earth-moon  barycenter,  refer- 
enced to  inertial  coordinate  system  (I,J,K) 

R - vector  from  earth-moon  barycenter  to  satellite, 
referenced  to  x,y,z-frame 

r - vector  from  earth-moon  barycenter  to  satellite, 
referenced  to  I,J,K-frame 

Substituting  the  above  expressions  into  equation  (2),  the 
equations  of  motion  can  be  written  as 

2 . ... 


_2  k^M  _ 

(U)x  + R + 2a)^R  + R = - ^ Tgp 

ep 


7 ” r + ■ A r ^ 

3 mp  m - sb 

mp  ^ 
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Defining  the  following  quantities 

F 

6F_  = 


and  (Ref  Figure  3) 


8 " Mp  ~ 

“n. 

V = - ® X 

r r 
ep  mp 


€ 

R = lij; 
P 


where 

6F  Solar  perturbative  force  on  satellite 
s 

V Gravitational  potential  fxmction 

M , - Mass  of  earth,  moon 

e m 

equation  (3)  can  be  expressed  as 

* • • • • • av  * 2 

€ - 2dl|;  - Gljj  = + 0 € + 6Fgg 

* * * • * * av  * ? 

* + 286  + 66  = ||  + 8“=*  + 

• • av 


(4) 


For  convenience  of  notation,  a potential  function  can 
be  written  as  (Ref  12 i 9) 

M M • n o o 

P = ^ + i 0^  (€^  + il;^) 

^ep  ^mp 

Using  this  notation  the  equations  of  motion  are 


• • • • 


dP 


€ - 2011;  - = if  + 6F 


s€ 


ap 


ij;  + 20€  + ©€  = + 6F^, 


sil; 


(5) 


0 = — + 6F 

° ao  “""so 


Equations  (5)  are  the  equations  of  motion  of  the  satel- 
lite with  respect  to  the  earth-moon  barycenter.  These  equa- 
tions include  the  following  assumptions.  First,  the  center 
of  mass  of  the  sun  is  the  center  of  mass  of  the  four  body 
solar  system.  Second,  the  earth  and  moon  are  assumed  to 
move  in  perfect  elliptical  orbits  about  their  barycenter j 
and  the  earth-moon  barycenter  moves  in  an  elliptical  orbit 
about  the  sun.  Third,  all  masses  are  assumed  to  be  point 
masses . 

One  can  also  note  that  equations  (5)  represent  the  super- 
position of  the  differential  solar  force,  6F  , and  the  three 

s 

body  equations  of  motion.  Thus,  the  derivation  can  be  com- 
pleted in  two  steps.  First,  the  three  body  equations  can  be 
linearized  about  the  point.  Next,  the  differential  solar 
force  can  be  linearized  about  L|^.  Once  these  expressions 
are  obtained,  they  can  be  substituted  into  equation  ($)  to 
give  the  linearized  four  body  equations  with  respect  to  the 
point. 

Three  Body  Equations  of  Motion 

The  three  body  equations  can  be  obtained  from  equation 

(5)  t>y  letting  the  differential  solar  force,  6F  , be  zero 

s 

(Ref  12 1 1) . This  gives 


€ - 29l|;  - SiJ/  = H 
+ 29€  + 9€  = 


The  vinlts  can  be  non-dlmenslonalized  by  choosing  the  mass 
unit  to  be  the  total  mass  of  the  earth  and  moon  (M^  ■*■**_  = !) 
and  the  distance  unit  to  be  the  earth-moon  distance  at  closest 
approach.  Taking  epoch  (t^  = 0)  to  be  the  time  of  closest 
approach,  Kepler's  law  gives 

k2(M  +3,)  , 


where 

N - mean  angular  rate,  2Tf/Period 

a - semi-major  axis  of  satellite  orbit 
2 

k - gravitational  constant 

2 

Choosing  the  time  scale  such  that  k and  the  lunar 

period  are  unity  results  in  a time  unit  equal  to  4.384 

2 T 


mean 


solar  days;  thus  N = — ^ . But  from  Fig,  4,  D = 1 = 


(j.-e)  and  thus  = fl  - e3^.  The  non-dimensional  equations 


can  then  be  written  as 


» • • • • • 1 

e - 264,  - 9*  = L 11 

» • • • • • ^ 

+ 29€  + 06  = ^ 

- i_  iP 
P - n2 


where 


0 is  the  true  anomaly  of  the  lunar  orbit 

t 

0 is  the  speed  of  the  rotating  frame  relative  to  the 
inertial  frame 

Using  Kepler's  law,  the  angular  velocity  of  the  rotating 
coordinate  system,  0,  can  be  reduced  to  (Ref  12:3-4) 


(1-e  Cos  (i)' 


Differentiating  gives  (Ref  12 i 3-4) 


when 


e is  the  lunar  eccentricity 
4 is  the  lunar  eccentric  anomaly 
For  small  perturbations  (x,y,z)  about  the  point 
(Xq*  Yq*  0),  the  coordinates  of  the  satellite  can  be  ex- 
pressed as  Xq  + X.  + yi  z (see  Figure  4).  The  right  hand 
side  of  equation  (7)  can  be  expanded  about  the  equilibrium 
point  (x  f y , 0)  in  a Taylor  series  as 


where 


Xq  = short  hand  notation  for  coordinates  of  the  lA 
point  (Xq,  y^,  o) 

It  is  assumed  that  the  perturbations  (x,y,z)  from  the 

point  will  be  small  compared  to  the  distance  from 

(Xo»  yQi  o)  to  the  earth-moon  barycenter.  Also,  the  terms 
2 2 2 

X , y , z , xy,  xz,  yz,  will  be  much  smaller  than  x,  y,  or  z 


and  can  be  neglected.  Only  the  zero  and  first  order  terms 
of  equation  (8)  will  be  retained. 

Making  use  of  the  following  relationships  (Ref  12 i 4-8) 


2 (1-e 


= (1  - e)^ 


where 

H is  the  eccentric  anomaly  of  the  lunar  orbit 
e is  the  eccentricity  of  the  lunar  orbit 
m is  the  mass  of  the  moon  (non-dimensionalized) 
the  partial  derivatives  of  equation  (8)  can  be  evaluated  at 
the  equilibrixim  point  (x^,  y^,  0)  as  (Ref  12*7  - 8) 


N~  de 


e Cos  n) 


e Cos  Li) 


e Cos  M. 


e Cos  u 


where 


5?! 


= 

— o 


After  evaluating  the  partial  derivatives  and  noting  that 
(Xq*  3^0’  a solution  to  equation  (7),  the  linearized 

equations  of  motion  about  the  point  are 


X - C 2^-e^  F^]  y + [2e  Sin  nVl-e^  F^]  y =C  (l-e^)F^  - iF^]^ 

+ QP3]  y 

y -+  C 2Vl-e^  F^3  * " C2e  Sin  u-A-e^  P^]  x = ^ 


Cd-e^jp'*  + 4 P3]  y 


(9) 


z = -F-'z 


where 


M 


Q = l-2itn  Mg  + = 1;  m = ; F = l/(l-e  Cos  \i) 


Restricted  Three  Body  Equations  of  Motion.  It  should 
be  noted  that  the  linearized  equations  of  motion  for  the  re- 
stricted three  body  problem  can  be  easily  obtained  from 
equation  (9)  by  letting  the  eccentricity,  e,  approach  zero. 
The  linearized  restricted  three  body  equations  are  (Ref  5: 
22-23) 


X - 2y 


= ^ X . ^ (1  - 2m)y 


- ^ X . 2 


• • 

z = -z 
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Four  Body  Eauations  of  Motion 


To  complete  the  derivation  of  the  four  body  equations 

of  motion  linearized  about  the  point,  an  expression  for 

the  differential  solar  force,  6F  , linearized  about  the  L,, 

S 

point  is  needed.  6F_  represents  the  difference  between  the 

s 

solar  gravitational  attraction  of  the  sxin  on  the  satellite 
and  the  solar  gravitational  attraction  on  the  earth-moon 
barycenter.  The  derivation  and  linearization  are  presented 
in  detail  in  Appendix  A.  The  linearized  expression  for  6F 

s 

is  summarized  here  as 


= -N 

X + X 

+ 3N 

0 

o 

+ 

z 

(Xq  +x) 


yo  + y 


(1  + Cos  2P)  - I Sin  23 


-(x  + x) 

g Sin  23  + g (1  - Cos  23 


-(yp  y) 
2 

2 

^yp  y> 

2 


^ ^o  yp 


(1  + Cos  20)  + 

Cos  3 (3N  Sin  20) 

Cos  3 + I Sin  3 
Sin  3 

^ S0  + ^ Cos  23  Sin  0 
Cos  3 Cos  0 
Sin  23  Sin  0 
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N 


\ 


e 

®] 

a 

Un 


= “/"sb 

= (1  - 2m)  (1  - e Cos  |i)/  (1  - e)  2 
= (1  - e Cos  n)/(2  (1  - e)) 

- lunar  orbit  eccentricity 

- earth-moon  barycenter  orbit  eccentricity 

- lunar  orbit  eccentric  anomaly 

- earth-moon  barycenter  orbit  eccentric  anomaly 
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Table  I 
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Coefficients  of  the  Equations  of  Motion 


Term 

Expression 

=1 

a/L-e^ 

<=2 

(l-e2)F^  -iF^  + 1 (3  cos 

2P  - 

1) 

•^JpQF^  -2e 

sin  4 Vl-e^  F^ 

“4 

^ sin  2e 

<=5 

f *0  (3  cos 

2P  - 1) 

06 

1 (1  + 

cos  2P} 

O7 

ny^  cos  P 
2 0 

(CP  - 5033)  - 

2^0 

CPl 

09 

='c’^0  (i 

+ ^ cos  2P) 

0 

H 

0 

^ N (1  + cos  2P) 

=11 

!>  n cos  P 

^12 

(cos  P - 5 cos 

33)} 

=13 

+ 5 cos  2P) 

^ -H—  X (1 
“ "sb  0 ‘ 

+ 5 cos  2P) 

^15 

.1 


I 

i 


i 


5^  *0  » 

^ QP^  + 2e  sin  s 

(1-6^)f'*  + J P^  + I 


2 5^0 


r^o 

|a  *0  oos  » 


2^  *0  ^o  “■>=  S 

5^  <3  *0  <1  * 5e  2S)  + yl) 


2^  ==0  --  ^ 

7-23—  y 
4 R3^  ^0 

¥t7  \ Si"  2S 

SD 


y.  Sin  2P 


f 


Term 


Table  I - continued 


Expression 


32 

-F^  -g  (1  + 3 cos 

33 

^ *0  Sin  23 

34 

Sin  23 

35 

^ Sin  3 

36 

(1  + cos  23) 

37 

38 

B;;  Vo 

39 

^16 

40 

y„  Sin  2B 

41 

8 *0  ® 

42  . 

5^yo 

9_L 

1 - 

e 

cos  u) 

2 

(1- 

e) 

a_ 

- e 

cos  u) 

2 

rr 

- e) 

1 - 

2m 

^b 

®i 

cos  a,  \ 

- 

®m^ 

1/(1  - e cos  4) 
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A further  simplification  of  the  equations  of  motion, 
equations  (11)  - (13),  can  be  made.  Looking  at  Table  I one 
can  see  a number  of  coefficients  where  the  factor 
appears  explicitly.  For  example,  typicsd  values  of  C2  and 
are  .0037  and  .000027,  respectively.  The  resultant 
equations  neglecting  these  terms  are 

X - (Cg  + Cos  2^f)^  - (C^  - Sin  20)^  - (-C^ 

Cos  20)^  = C^  + C^  Cos  2j2f  - C^  Sin  20  (14) 

y’  + Cj^  * - iC^Q  - Sin  20)^  - (C^^  - C23  Cos  20)^ 

- Sin  20)^  = C20  - C21  Cos  20  - Sin  20  (15) 

V - (-Cj^  -Cj^  Cos  20)^  - (C^^  Sin  20)^  - (0^2'*’  ^36 
= -C^^  -C^^  Cos  20  + C^j^  Sin  2)2(  (I6) 

where  the  coefficients  are  as  listed  in  Table  I. 

Very  Restricted  Four  Body  Problem 

Just  as  the  linearized  equations  of  motion  for  the  re- 
stricted three  body  case  could  be  obtained  from  the  three 
body  case,  the  very  restricted  four  body  equations  can  be  ob- 
tained from  the  equations  of  motion  for  the  four  body  case 
(Eqn  14-16) . The  very  restricted  four  body  problem  implies 
that  the  orbits  of  the  earth  and  moon  about  their  barycenter 
and  the  earth-moon  barycenter  about  the  sun  are  circ\ilar 
(i.e.,  e=  ej^=  0).  Another  restriction  is  that  all  orbits  are 
coplaner  (i.e.,  3=0).  When  these  conditions  of  the  very 
restricted  four  body  problem  are  imposed,  the  equations  re- 
duce to  those  developed  by  deVries  (Ref  4:18-19)  who  started 
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with  the  very  restricted  four  body  model.  Figures  5 and  6 
provide  a comparison  between  the  trajectory  and  control 
determined  by  the  very  restricted  four  body  equations  (Ref  8) 
and  the  trajectory  and  control  as  determined  by  the  four  body 
equations  (Eqn  14-16)  with  the  limar  and  barycenter  eccentric- 
ity zero  (e  = e^^  = 0)  and  the  inclination  of  the  ecliptic 
zero  O = 0).  The  similarity  of  the  two  verifies  the  four 
body  derivation. 

Assumptions 

The  following  is  a summarization  of  the  assumptions  used 
in  the  derivation  of  the  linearized  equations  of  motion  about 
the  point  for  the  four  body  system  of  Figure  2. 

1.  All  masses  are  point  masses.  The  center  of  mass  of 
the  sun  is  assumed  to  be  the  center  of  mass  of  the 
four  body  solar  system. 

2.  The  orbits  of  the  earth,  moon  and  earth-moon  bary- 
center are  assumed  to  be  elliptical. 

3.  The  fourth  body  (satellite)  is  small  enough  that  it 
does  not  affect  the  motion  of  the  other  three  massive 
bodies. 

4.  The  perturbation  of  the  satellite  from  is  small 
compared  to  the  distance  from  to  the  earth-moon 
barycenter. 

Normalized  Units 

Unless  otherwise  noted  all  quantities  are  expressed  in 
normalized  units.  The  astronomical  constants  upon  which  the 
normalizations  are  based  are 
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< Figure  5 

I Trajectory  - Very  Restricted  Four  Body  Equations  (Ref  808) 


I 


Figure  6 

Trajectory  - Four  Body  Equations  (e=  ej_  = B = 0) 


PTT 


1/m  = 82.301 
m„  = 328905.2 
n = .0748 

(I)  =1 

= 1 

The  conversion  from  normalized  units  to  physical  units 
can  be  made  through  the  following  relationships  (Ref  li429) 

1 time  unit  = 4.384  mean  solar  days 

1 distance  unit  = 384400  kilometers 

1 velocity  unit  = 1023  meters/second 

2 

1 acceleration  unit  = .002724  meters/second 
Known  Station  Keeping  Costs 

In  the  previous  works  where  the  plauier  problem  has  been 
considered,  the  cost  of  station  keeping  has  vauried  from  zero 
(for  the  vincontrolled  satellite  (Ref  11:84))  to  .OO8567  nor- 
malized acceleration  units  (for  the  perfectly  controlled 
satellite  (Ref  8:6-7,  29-31))- 

The  cost  for  perfect  station  keeping  can  be  fovind  by 
using  a control  which  exactly  cancels  the  forcing  terms  of 
the  equations  of  motion,  (Equations  (14)  - (I6))  i.e., 


= -(C^  + cos  2/  - Cr,  sin  20) 

Uy  = -(CgQ  - cos  20  - C22  sin  20) 

20  + sin  20) 

The  average  magnitude  of  the  thrust  acceleration  is 


(17) 


= (1/2tt)  + uj  + U^)  d0 


(18) 
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But  for  the  case  at  hand  the  coefficients  of  the  forcing 
terms  of  equation  (18)  are  time  varying.  Considering  a 
constant  co  and  Q as  in  the  very  restricted  four  body  case 
the  variable  of  integration  can  be  transformed  to  dt  by 
the  relations 

= a + (u)  - Q)t  j d^(  = (u)  - n)  dt 
Thus,  equation  (18)  becomes 

“av  ‘“x  * "y  * 

Using  the  computer  to  integrate  the  equation,  the 
average  thrust  acceleration  for  zero  drift  is  .008697  norm- 
alized acceleration  units.  (Note  that  the  limits  of  0 from 
O-Ztt  correspond  to  time  limits  of  0-6.82  normalized  time 
units) . 


Performance  Index 

It  is  useful  to  define  a quauitity  that  will  given  an 
indication  of  the  cost  of  controlling  the  satellite  in- 
dependent of  the  vehicle  weight  or  rocket  performance . 
Looking  at  the  rocket  equation 

[Fuel  weight]  = [Final  weight]  [exp  (^V/(g  I )-l)]  (20) 

o sp 

one  can  see  that  minimizing  the  total  velocity  increment, 

AV,  will  minimize  the  fuel  required  for  a given  specific 

impulse  (I„_)  and  no  staging.  So,  AV  seems  to  be  an  apt 
sp 

quantity  to  use  as  a performance  index.  The  total  velocity 
increment,  AV,  can  be  catLculated  by 


(21) 


n ^ 


AV  = Xj  |u(t)|  dt 

where  U(t)  is  the  control  specified  by  the  quadratic 
cost  function,  equation  (28). 
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II.  Analysis 


Optimal  Control  Problem 


For  a linear  system,  the  optimal  control  with  respect 
to  a quadratic  cost  function  can  be  obtained  as  a linear 
combination  of  the  states  (Ref  3«  148-151). 

Given  the  linear,  time  varying  system 


X(t)  = A(t)  X(t)  + B(t)  U(t) 


and  the  quadratic  cost 


J = i X(T)^  X(T)  + i/g  (X(t)’^  Q(t)  X(t)  + UCt)"^  R(t)  U(t))  dt 


where  and  Q(t)  are  positive  semidefinite  matrices  and 
R(t)  is  a positive  definite  matrix. 

The  optimail  control  is  given  by 

U(t)  = -R(t)’^B(t)'^K(t)  X(t)  (24) 

where  K(t)  is  the  unique  solution  to  the  matrix  Riccati 
equation 

K(t)  = -K(t)  A(t)  - A(t)'^  K(t)  + K(t)  B(t)  R(t)"^  B(t)'^  K(t) 

-Q(t)  (25) 

subject  to  the  boundary  condition  K(T)  = 

Figure  7 gives  the  form  of  the  optimal  lystem.  The  state 
equations  (Eq  21)  can  be  written  as 

X(t)  = [A(t)  -B(t)  R(t)'^  B(t)'^  K(t)]  X (t)  (26) 

For  this  system  the  matrices  R(t)  and  B(t)  are  specified 
as  constants.  Thus,  the  matrix  K(t)  determines  the  gain  of 
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of  the  feedback  path  and  the  behavior  of  the  system.  K(t) 
will  be  referred  to  as  the  "gain"  matrix.  As  seen  from 
equation  (25) • it  is  independent  of  the  states.  Once  the 
matrices  A(t),  B(t),  R(t),  Q(t),  and  the  final  time  are 
specified  equation  (25)  is  solved  backwards  in  time  to  give 
the  time  varying  gain  schedule. 

State  Equations 

The  four  body  equations  of  motion  for  the  satellite  with 
respect  to  the  L4  point  were  derived  in  chapter  I and  are 
restated  here  for  convenience  as 

• • • 

X = y + (Cg  + cos  20)  X + (C^  - sin  20)y 

+ Cos  20)  z + Cos  20  - Cg  Sin  20  (l4) 

y = i + (C^  -C^  Sin  2^)x  + (C^^q  Cos  2^)y  + (C^^  Sin  20)  z 

+Ci2  - Cos  20  - C^^  Sin  20  (15) 

z = (“Cj^^  Cos  20) X + (Cj^g  Sin  20)y  + (CgQ  + Cos  20)  z. 

“^17  ”^17  ^ ^18 

where  the  coefficients  have  been  reniunbered  and  aire  listed 
in  Table  II.  The  equations  of  motion  can  be  put  in  the  form 
of  equation  (25)  choosing  the  following  states 

» • • 

x^  = x;  x^  = xj  x^  = y;  x^  = yi  x^  = 1;  x^  = z;  x^  = zj 

Xg  = cos  20 t x^  = sin  20 
The  state  equations  are 
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i 
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Table  II 

Coefficients  of  the  State  Equations 


Term 

Expression 

f 

2>/l-e^ 

1 

=2 

(l-e^)F^  - + 1 (3  cos  23  - 1) 

=3 

^ N (1  + cos  23) 

^ QF^  - 2e  Sin  F^ 

=5 

^ N cos  3 

C6 

J Xq  (3  cos  23  - 1) 

=7 

^ N x^  (1  + cos  23) 

2 N Yq  cos  3 

r 

1 

^ QF^  + 2e  sin  n7l-e^  F^ 

i 

f 

i 

r 

=10 

(l-e^jP^  + ^P3  * f 

i 

i 

I * 1 

‘'ll 

3N 

2 

i ■ 

[ ' 

Cl2 

^ y 

2 ^0 

t 

i 

i , 

=13 

^ y 

2 ^0 

( 

\ 

‘'l^ 

^ x„  cos  3 

2 0 
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*1  *2 

Xg  = (Cg  + C3  cos  20)  + (C^j^  -C3  sin  2^)x3  + 

+ cos  20)  x^  + C^Xg  - CgX^ 

X3  = 

Xji^  = (C9  - C3  Sin  20)  Xj^  "^2*2  '*’  ^^10  " ^11  °®®  2^)  *3 

■*■  ®12  *5  ^19  *6  " ^13^8  - ^14*9 

X3  = 0 (27) 

^6  = *7 

X^  = -C^3  cos  20)  x^  + C^g  Sin  20  x^  - x^ 

'*’  ^^20  *^21  °°®  ^^^*6  " ^17*8  '*’  ^12*9 

Xg  = -20  X9 
X9  — 2o  Xg 

The  object  is  to  minimize  drift  and  fuel  used;  thus,  the 
unconstrained  quadratic  cost  function  used  is 

J = i /q  Cq  (x^  + x^  + x^)  + (U^  + u|  + u|)]  dt  (28) 

This  is  in  the  form  of  equation  (23),  where 


K-  = [0]  (a  9x9  zero  matrix) 


M2- 


1 


Using  the  matrices  given  above  and  performing  the  In- 
dicated operations,  Eqn  (24)  and  (22)  can  be  simplified  to 


^21  ^22  ^23  ^24  ^25  ^26  ^27  ^28  ^29 


U(t)  — - ^4l  ^42  ^43  ^44  ^45  ^46  ^47  ^48  ^49  — (29) 

^71  ^72  ^73  ^74  ^75  ^76  ^77  ^78  ^79 


10  0 0 0 0 0 0 0 0 


^21  ^22  ^23  ^24  ^25  ^26  ^27  ^28  ^29 


0 00000  0 00 


^41  *^42  ^43  ^44  ^45  ^46  ^47  *^48  *^49 

x(t)  = A(t)-  000000000  x(t)  (30) 


0 00000000 


^71  ^72  ^73  ^74  ^75  ^76  ^77  ^78  ^79 


0 00000000 


0 00000000 
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III.  Gain  Matrix  Computations 


The  gain  matrix,  K(t),  is  computed  by  solving  the  matrix 
Riccati  equation  (Eq  25)  backwards  from  the  final  time  (T). 
The  solution  is  optimal  gain  schedule  based  upon  the  spec- 
ified final  time.  This  gain  schedule  would  then  be  used 
to  control  the  satellite  as  specified  by  the  control  equation 


U(t)  = -R(t)'^  B(t)’^  K(t)  X(t) 


Figures  8-10  present  plots  of  typical  gain  elements 
versus  time  for  a variety  of  weights,  q,  (Ref  Equation  (28)). 
As  seen  in  the  plots  all  exhibit  common  factors.  Each  has 
a trzmsitory  phase  and  a steady-state  phase  that  is  periodic 
with  approximately  the  lunar  period.  As  the  weighting,  q, 
is  increased,  the  transitory  phase  shortens,  and  the  peak 
amplitude  in  the  transient  phase  decreases.  The  period  re- 
mains unchanged  for  the  particular  gain  element. 

In  the  matrix  of  coefficients  of  the  states,  the  "A" 
matrix,  there  is  a dependance  upon  the  eccentric  anomalies 
of  the  earth-moon  barycenter  and  the  moon.  This  leads  to 
the  question  of  what  effect  the  relative  positions  of  the 
earth,  moon  and  sun  have  upon  the  computations.  Figures  11- 
13  depict  gain  elements  computed  for  q = 1 for  a variety  of 
initial  conditions.  The  changes  from  one  plot  to  another  are 
quite  small,  but  the  effects  on  system  performance  are  signif- 
icant. Chapter  IV  will  examine  the  effects  of  initial  con- 
ditioning in  detail. 
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Figure  11 


Gain  Element  Vs  Time 
(q  = 1;  a = 30°) 
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Figure  12 

Gain  Element  Vs  Time 


GAIN 


Constant  Gain  Elements 


While  the  varying  gain  schedule  above  gives  the  optimvun 
path  and  control  for  the  satellite,  it  is  a very  costly  method 
in  terms  of  computation  and  storage  requirements  and  so,  there  j 

is  motivation  to  derive  a less  cumbersome  gain  schedule.  Since  j 

the  gains  are  periodic  in  the  steady-state,  averaging  the  I 

i 

gain  over  the  steady-state  interval  should  provide  a near  ' 

optimal  constant  gain.  The  use  of  a constant  gain  has  several  | 

j 

advantages.  First,  there  is  no  requirement  for  a complete  ? 

i 

solution  to  the  Riccati  equation.  It  is  only  necessary  for  j 

the  Riccati  equation  to  be  solved  for  a short  period  so  that 
the  constant  gain  can  be  computed}  thus,  there  is  no  depend- 
ence on  final  time  as  there  is  for  the  optimal  gain.  Second, 
the  only  storage  required  is  for  a single  gain  matrix  rather 
than  an  entire  gain  history.  It  must  be  noted,  however,  that 
the  price  for  the  simplicity  of  the  constant  gain  is  sub- 
optimal  station  keeping.  Table  III  lists  the  constant  gains 
computed  for  a variety  of  weights  with  a = 0.  As  the  weight- 
iriSi  <l»  increases  the  values  of  the  gain  increase.  In  terms 
of  performance  it  will  be  shown  in  Chapter  IV  that  this  in- 
crease in  gain  values  will  result  in  a smaller  variance  from 
the  L4  position  and  an  increase  in  cost.  Table  IV  shows  the 
variance  of  the  gains  as  the  initial  sun  direction,  a,  is 
varied . 

The  computer  programs  used  to  compute  the  fixed  and  vary- 
ing gains  are  listed  in  Appendix  B. 

; 
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Table  III 


Fixed  Gains 


.01 

-6.6937E-01 

.02 

-7.3OO3E-OI 

.03 

-7.4116E-01 

.10 

-7.0746E-01 

.30 

-4.9936E-OI 

1.00 

3.2440E-02 

3.00 

9.2753E-OI 

10.00 

2.5732E+OO 

100.00 

9 . 6678E+OO 

l.E+04 

9.9749E+OI 

5.5528E-OI 

6.0598E-OI 

6.319OE-OI 

7.IO2IE-OI 

8.2081E-01 

I.O738E+OO 

1.4897E+00 

2.19I8E+OO 

4.3IO5E+OO 

1.4112E+01 


I.I335E+OO 
1.2463E+00 
-1.2954E+00 
-1 . 3585E+OO 
-I.3I73E+OO 
-I.266IE+OO 
-1 . 2288E+00 
-I.5I7IE+OO 
-3.2773E+OO 
-1.2855E+01 


.01 

8.2675E-01 

-4.5763E-01 

.02 

I.O583E+OO 

-5.I332E-OI 

.03 

1.2093E+00 

-5.4028E-01 

.10 

1 . 725IE+OO 

-5.4847E-01 

.30 

2.2759E+00 

-4.6732E-OI 

1.00 

2.93IOE+OO 

-2.8703E-OI 

3.00 

3.5645E+OO 

-8 . 7659E-02 

10.00 

4.2697E+00 

8.6078E-02 

100.00 

6.0692E+OO 

I.7286E-OI 

l.E+04 

1.5434E+01 

8.1488E-02 

I.2747E+OO 
I.5912E+OO 
1.7877E+00 
2.3828E+OO 
2.9463E+OO 
3.558IE+OO 
4.1924E+00 
5.3132E+OO 
I.I592E+OI 
I.OI3OE+02 


q 

K^l 

.01 

2.0049E-04 

.02 

I.639OE-O4 

.03 

1.2442E-04 

.10 

-4.7705E-05 

.30 

-2.1486E-04 

1.00 

-3.8574E-04 

3.00 

-4.9425E-04 

10.00 

-5.5786E-04 

100.00 

-5.9674E-04 

l.E+04 

-6.3627E-04 

-2.3693E-O 
-2.33I8E-04 
-2.2747E-04 

-2.1019E-04 

-2.0717E-04 

-2.1272E-04 

-2.0856E-04 

-1.8841E-04 

-I.2623E-04 

■4.4773E-05 


.1878E-0 


3. 9 70 2E -04 
3.6286E-04 
2.I3OOE-04 
8.1473E-05 
2.2487E-06 

I.2445E-05 

I.7092E-06 

2.7829E-05 

4.8417E-05 


Table  III 


continued 


q 

''a-* 

''25 

NO 

CM 

.01 

-4.5763E-OI 

-I.I067E-03 

-2.342IE-O4 

.02 

-5.I332E-OI 

-I.2IO9E-03 

-2.3417E-04 

.03 

-5.^028E-01 

-I.2509E-03 

-2.329IE-04 

.10 

-5.4847E-01 

-I.2619E-03 

-2.2460E-04 

.30 

-4.6732E-OI 

-1.093^-03 

-2.3102E-04 

1.00 

-2.8703E-OI 

-6.5882E-04 

-2.6752E-04 

3.00 

-8.7659E-02 

-7.3255E-0^ 

-3.4941E-04 

10.00 

8.6078E-02 

6.0977E-04 

-4.4524E-04 

100.00 

I.7786E-OI 

I.22I8E-O3 

-5.5599E-04 

I.E+04 

8.1488E-02 

I.23OIE-03 

-6.3265E-04 

^45  ^46 


I.2S80E-03 

2 . 7118E-04 

.02 

9.0443E-01 

I.5905E-03 

3.063IE-04 

.03 

1.0466E+00 

I.7927E-03 

3.3044E-04 

.10 

I.4969E+OO 

2.433IE-03 

3.0784E-04 

.30 

1.9446E+00 

3.03IIE-03 

2.0121E-04 

1.00 

2.4173E+00 

3.6OO8E-03 

5.5842E-05 

3.00 

2.8255E+OO 

3.8814E-03 

-4.6689E-05 

10.00 

3.3068E+OO 

3.7539E-03 

-9.5309E-05 

100.00 

4.8853E+OO 

2.8743E-03 

-7.4269E-05 

l.E+04 

1.4246E+01 

2.2575E-03 

-I.O638E-05 

^74  ^75  ^76 


.01 

8.2222E-05 

-5.7368E-O6 

3.8978E-03 

.02 

7.3^97E-05 

-7.828OE-O6 

8.4148E-03 

.03 

5-5020E-05 

-9.433OE-O6 

I.3322E-02 

.10 

-3.2374E-O5 

-I.8953E-O5 

4.5687E-02 

.30 

-I.0709E-04 

-4.26IOE-O5 

I.3455E-OI 

1.00 

-I.393IE-04 

-9.6705E-05 

4.0404E-01 

3.00 

-1.2473E-04 

-I.635OE-04 

9.908OE-OI 

10.00 

-8.875IE-05 

-2.2822E-04 

2.3027E+OO 

100.00 

-3.2923E-05 

-2.9^53E-04 

9.O302E+OO 

l.E+04 

-I.6929E-O6 

-3.27IOE-04 

9.8997E+OI 

Table  III 


continued 


q 

^27 

IN} 

CD 

*^29 

.01 

-2.3693E-o4 

-6.2699E-04 

-3.7076E-03 

.02 

-2.3318E-04 

-4.9653E-04 

-4.232OE-O3 

.03 

-2.2747E-04 

-4.3649E-04 

-4.5232E-O3 

-5.4612E-03 

.10 

-2.1019E-04 

-3.3435E-04 

.30 

-2.0717E-04 

-1.7917E-04 

-6 . 7004E-03 

1.00 

-2.1272E-04 

7.650IE-04 

-8.I93OE-O3 

3.00 

-2.0856E-04 

2.4143E-03 

-8.6067E-03 

10.00 

-1.8841E-04 

3.2614E-03 

-7.8685E-O3 

100.00 

-1.2623E-04 

3 . 643IE-03 

-6.5406E-03 

1 .E+04 

-4.4773E-05 

3.6056E-03 

-6.5789E-O3 

q 

K47 

^48 

K49 

.01 

8.2222E-05 

-I.7533E-O3 

4.I33IE-O3 

.02 

7.3497E-O5 

-2.9406E-03 

4.8988E-03 

.03 

5.502OE-O5 

-3.815IE-03 

5.2554E-O3 

.10 

-3.2374E-05 

-7.O729E-O3 

5.4280E-03 

.30 

-1.0709E-04 

-1.0074E-02 

3.8353E-O3 

1.00 

-I.393IE-04 

-I.I522E-O2 

5.8308E-04 

3.00 

-1.2473E-04 

-1.0809E-02 

-I.7835E-O3 

10.00 

-8.875IE-05 

-9.6I66E-03 

-2.8655E-O3 

100.00 

-3.2923E-05 

-8.1452E-03 

-3.8II6E-O3 

l.E+04 

-I.6929E-O6 

-6.7OI2E-O3 

-3.8I87E-O3 

Table  IV 
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IV.  System  Performance 

Optimum  versus  Modified  Control 

To  make  the  station  keeping  costs  of  value  to  the  user, 
system  performance  must  be  evaluated  in  terms  of  cost  versus 
distance  from  the  lA  point.  In  chapter  III  both  constant  and 
time  varying  gains  were  computed.  These  lead  to  two  types 
of  controllers!  the  optimum  controller,  using  the  varying 
gain  schedule;  and  the  modified  controller,  using  the  constant 
g£dn.  Table  V gives  a comparison  of  the  optimum  versus  mod- 
ified controller  for  q = 1.  (Note:  Quadratic  and  AV  costs 
are  computed  by  equations  (20)  and  (21)  respectively) 


Table  V 

System  Comparison 


Final  Time 

System 

Peak  Drift 

AV  Cost 

Quadratic  Cost 

15 

Optimal 

.0046075 

.064284 

.00093150 

15 

Modified 

.0050679 

.066905 

.00092397 

25 

Optimad 

.0046075 

.103100 

.0015483 

25 

Modified 

.0050679 

.10589 

.0015445 

45 

Optimal 

.0046075 

.17738 

.0026726 

45 

Modified 

.0050679 

.17635 

.0027804 

If  AV  cost  were  the  only  consideration,  the  optimal  con- 
troller would  be  the  choice.  However,  when  one  considers 


computational  requirements,  the  modified  controller  is  better 
in  that  it  provides  near  optimal  performance  with  a fraction 
of  the  computational  requirements  of  the  optimal  controller. 
Table  V shows  that  as  time  increases  the  actual  costs  of  the 
two  systems  become  nearly  equsd  while  the  difference  in  peak 
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drift  Is  constant.  Figures  14  and  15  present  graphically 
the  satellite  path  versus  time  for  the  optimal  and  modified 
control  systems  respectively. 

One  can  note  from  the  figures  that  both  controllers 
exhibit  a transient  period  starting  at  x=y=z=0.  The  optimal 
controller  drives  the  satellite  to  the  steady  state  in  approx- 
imately 2.5  time  xmits;  the  modified  controller  takes  approx- 
imately 5 time  units  to  reach  steady  state. 

Constant  Gain  Performance 

For  the  remainder  of  this  study,  the  constant  gain 
(modified)  controller  will  be  examined.  The  results  will 
hold  for  the  optimal  controller  as  well.  Appendix  C lists 
the  computer  program  used  for  either  constant  gain  or  optimal 
gain  schedules. 

In  chapter  II  it  was  noted  that  there  were  two  extremes 
of  costt  perfect  station  keeping  and  no  control.  But  what 
of  the  region  between  the  two  extremes?  Figure  I6  depicts 
the  path  of  the  satellite  for  an  approximate  one-year  mission 
when  no  control  is  used.  The  resxiltant  path  is  unstable. 
Figures  1?  - 19  depict  the  path  for  aji  approximate  one-year 
mission  when  various  weights,  q,  are  used.  In  all  cases, 
after  a brief  transient  period  the  satellite  enters  roughly 
a steady-state  path  and  remains  locked  in  that  path  as  long 
as  the  controller  is  operational.  Several  characteristics 
can  be  noted  from  the  figures.  First,  as  q is  increased  the 
duration  of  the  transient  phase  is  shortened.  Second,  the 
steady-state  path  becomes  tighter  with  increasing  q.  Fig- 
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Figure  14 

Satellite  Trajectory  - Optimal  Control 
(q  = 1;  a = OO) 
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SATELLITE  TRAJECTORY 


Figure  18 

Satellite  Trajectory 
(q  = 1.0;  zero  IC) 

1 

I 

I I 

,! 
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Figure  19 

Satellite  Trajectory 
(q  = 10. ; zero  IC) 
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1 

ure  20  gives  a typical  plot  of  the  velocity.  It  exhibits 

the  same  characteristics  as  the  trajectory.  As  can  be  noted 

from  equation  (24)  the  control  is  formed  by  a combination 

of  states,  thus  in  the  steady-state  the  control  is  periodic 

with  a nearly  constant  magnitude.  ] 

I 

The  performance  of  the  modified  control  can  be  summarized 

i 

by  Figure  21.  The  cost  vauries  almost  linearly  from  perfect  i 

control  down  to  q = .1  and  then  asymptotically  approaches  * 

zero  as  q is  decreased.  While  perfect  station  keeping  is  the 
most  costly  option,  controlling  with  q less  than  .1  is  not 
desired  because  of  the  large  increase  in  peak  radius  with 
minimal  reduction  in  cost.  Figure  21  will  enable  the  user 
to  choose  the  best  option  for  his  situation  by  evaluating 
cost  versus  tolerable  distance  from  L4. 

Effects  of  Initial  Conditions 

In  chapter  III  it  was  shown  that  the  gains  vary  with 
initial  configuration  of  the  earth,  moon  and  sun.  In  this  ^ 

section, • system  performance  will  be  studied  in  response  to 
initial  planet  configurations  as  well  as  initial  position 
and  velocity  conditions  for  the  satellite.  Figure  22  depicts 
the  range  of  average  thrust  as  a is  varied  while  initial 
satellite  position  and  velocity  are  zero.  One  can  see  that 
there  are  initial  conditions  that  will  minimize  the  cost  for 
a given  q.  Figures  23  - 24  give  an  indication  of  why  the  cost 
varies.  For  the  smaller  a the  transient  interval  is  longer,  ; 

i.e.,  the  controller  is  scheduling  a smaller  thrust  for  a | 

longer  period  of  time.  Thus,  the  cost  is  smaller.  Note  I 


58 


-.DO  HI.H9  92.3*1  i2H. 

NRUTICHL  MILES  )<ia 


MRX  IISTBr^JEE  rR3M  L4- 


Figure  21 


Total  Cost  VS  Peak  Drift 
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that  for  a given  q,  the  steady-state  path  remains  unchanged 
as  a is  varied.  By  proper  selection  of  ini tied  position 
and  velocity,  the  transient  phase  can  be  eliminated)  but, 
the  total  cost  is  higher.  This  is  because  the  controller  is 
initially  scheduling  higher  thrust  levels  than  it  is  in  the 
corresponding  transient  time  period.  From  Figure  22  one  can 
alsc  observe  that  as  q is  increased  the  dependence  on  initial 
conditions  is  decreased.  This  is  a direct  result  of  the 
shorter  transient  interval  associated  with  the  higher  values 
of  q.  Note,  that  for  q = 1,  average  thrust  varies  9%  and 
for  q = 10  thrust  varies  2.8%  (Ref  Figure  22). 
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V.  Discussion 


Model  Deficiencies 

The  three-dimensional  model  used  for  the  four  body 
system  is  restricted  by  the  following  assumptions! 

1.  Only  the  gravitational  perturbations  from  the  sun, 
earth  and  moon  were  considered.  Perturbative  effects 
such  as  solar  wind,  gravitational  pull  from  other 
planets,  etc.  were  neglected. 

2.  All  masses  were  considered  point  masses.  The  earth 
and  avn  were  considered  to  be  the  centers  of  mass 
for  their  respective  system. 

3«  The  nonlineeir  equations  were  linearized  about  the 
L4  point . 

Since  the  controller  used  is  formed  by  a combination  of  the 
states,  it  is  felt  that  model  deficiencies  were  minor. 

System  Realizability 

Thrust  Requirements.  The  thrust  requirements  for  station 
keeping  are  quite  small.  As  an  example,  (Ref  Fig.  21),  a 
1000-lb  satellite  would  require  .0017-lb  thrust  to  confine 
its  peak  drift  to  a 500  nautical  mile  radius.  So,  for 
continuous,  low  thrust  requirements  an  ion  or  monopropellant 
engine  would  be  the  most  likely  choice.  Since  the  thrust 
vector  is  periodic,  higher  thrust  engines  which  would  fire 
a si)ecified  discrete  times  might  provide  a thrust  which  would 
approximate  the  continuous  thrust  provided  by  the  modified 
controller. 
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Sensor  Requirements.  The  control  useu  requires  con- 
tinuous knowledge  of  the  position  and  velocity  states  and 
the  sun  direction,  0.  Laser  ranging  might  be  one  method  of 
providing  the  control  system  with  this  information.  The 
sensor  requirements  should  not  be  a serious  problem.  The 
constant  feedback  gains  would  have  to  be  calculated  a 
priori  based  on  the  desired  launch  date  (defining  initial 
solar  system  configuration)  and  the  best  estimate  of  the 
astronomical  constants  over  the  mission  duration. 


VI.  Si 


id  Recommendations 


The  linearized  equations  of  motion  for  a satellite 
stationed  at  or  near  the  lA  lunar  libration  point  were  de- 
rived for  the  three-dimensional I four  body  case.  A control 
system  was  then  obtained  by  employing  optimal  control  theory 
to  this  system  of  equations.  Modification  of  the  optimal 
controller  using  fixed  gains  resulted  in  a control  system 
that  was  more  computationally  attractive  while  providing  per- 
formance near  the  optimal. 

Several  observations  can  be  made  from  the  three-dimen- 
sional study  of  the  controlled  satellite.  First,  both 
trajectory  and  control  become  periodic  after  a brief  transi- 
tory phase.  This  allows  the  definition  of  average  thrust 
as  the  average  magnitude  of  the  control  over  one  period. 
Second,  initial  conditions  were  shown  to  affect  only  the 
transient  performance.  Third,  as  q is  increased  the  effects 
of  initial  solar  configuration  are  reduced,  i.e.  smaller  peak 
drift  but  longer  thrust.  Fourth,  initial  conditions  which 
eliminate  the  transient  phase  result  in  a more  costly  tra- 
jectory than  the  trajectory  with  a transient  phase.  And 
finally,  the  range  of  station  keeping  costs  is  bounded.  There 
is  both  an  upper  limit  and  a lower  practical  limit.  This 
enab3.es  the  user  to  make  the  tradeoff  between  maximum  thrust 

I 

he  can  afford  and  peak  drift  from  the  L4  point. 


I 


To  finally  validate  the  results  of  this  study  and  pre- 
vious ones*  the  control  system  should  be  tested  using  the 


nonlinear  equations  of  motion.  Using  actual  ephemeris  data 
and  comparing  the  results  with  this  and  previous  studies  will 
given  an  indication  of  the  accuracy  required  for  the  astro- 
nomical constants. 

Another  area  reqviiring  further  study  is  the  effects 
of  initial  conditions.  It  would  be  useftil  to  determine  a 
method  to  predict  the  optimal  initial  conditions  for  a given 
weighting,  q. 

Additionally,  a sensitivity  analysis  to  determine  the 
accuracy  required  in  determination  of  the  states  by  the 
sensors  is  necessary. 


E 
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Appendix  A 


Solar  Perturbative  Force 


The  solar  perturbative  force  is  given  by  (See  glossary 
page  vii  for  list  of  symbols  used) 


but 


_ P 1' 

6P  = — S-  - r 
®'s  mp  ^sb 


• * F. 

- b _ s - 

^sb  ■ M „ 3 ^sb 


(31) 


M r 


sb 


where  m is  the  mass  of  the  sxm  and  M is  the  mass  of  the  earth 
s 

and  moon.  Substituting  this  expression  into  equation  (31) • 
the  force  becomes 

J _ 

(32) 


^ ^ ;f  ^sp  - -sb3 


‘sb 


sp 


It  can  be  shown  that  (Ref  7:84-89) 

3 2 r„,  - r 2 -3/2 

= (1  + — — + V) 


(33) 


sp 


sb 


sb 


Using  the  expansion 

(1  + z)'"  = 1 - nz  + - nLn.±  W.Ln.i.Ap  + ... 


the  right  side  of  equation  (33)  becomes 


-5  3 . r 

( 3b\ 3 _ i sb 

^r__‘'  ~ ^ _2 


sp 


sb 


3 r-2  ,13  ‘■'sb  • 

2 2 U 

^ ^sb  ^sb 
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Equation  (3I)  then  becomes 

f?  .a-  - 3£ 

O l-^Sb  2 


3^sb‘^  3r^  , 

^8b  <l>  2 4 


8 7" 

^sb 


■‘'  * • ^’^sb  ■‘’  (35) 

Equation  (35)  is  a vector  equation  with  all  vectors 
taken  with  respect  *0  the  rotating  x,y,z-frame  (see  Figure  2). 
To  transform  the  vectors  from  the  I,J,K-frame  to  the  rotating 
frame  the  following  relations  are  used: 


Cos  0 Cos  (0  + ^)  Cos  0 Sin  (9  + 0) 

A = - Sin  {9+0)  + Cos  {9+0) 

- Sin  0 Cos  {9  + 0)  - Sin  P Sin  {9  + 0) 


Sin  P 
0 

Cos  P 


1 —X  _ aX  pi  _ p 

^sb  " ^sb  " ^sb 

Cos  P Cos  0 

1 r = 

€ 

- Sin  0 

1 

- Sin  P Cos  0 

- 

P 

where 

0 is  the  initial  sun  direction  relative  to  the  x,yjZ- 
frame 

P is  the  inclination  of  the  ecliptic  plane 
To  evaluate  the  vector  relationships  of  Equation  (35) 
the  following  are  useful 

I 

r ^ . r = R . (€  cos  P cos  0 - sin  0 - p sin  p cos  0) 
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• r)  =s  Rg^  (€  COB  3 cos  0 ~ sin  0 '■  p sin  3 cos  0) 

cos  3 cos  0 
-sin  0 
-sin  3 cos  ^ 

r)  r = Rg^  (€  cos  3 cos  0 - \jf  sin  0 - p sin  3 cos  0)  • 


Using  the  above  expressions,  Eqn  (35)  is 

N - € + 3(^  cos  3 cos  0 - sin  0 - p sin  3 cos  0)  • 

^ cos  3 cos  0 

P - sin  0 

- sin  3 cos  0 


cos  3 cos 


sin  3 cos  0)  ^ 
>1^ 


cos  3 cos  0 

- sin  0 

- sin  3 cos  0 


cos  S cos 


sin  3 cos 


cos  3 cos  0 I 

- sin  0 

- sin  3 cos  0 
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Appendix  B 


O FORTRAN  Program  for  Gain  Computations 

The  program  listed  in  this  appendix  is  used  for  the 
time  varying  gain  computations.  To  use  the  program,  the 
user  must  attach  and  library  the  IMSL  package,  CC66OO. 

This  provides  access  to  the  integration  routine,  ODE. 

The  output  is  both  printed  and  on  tape.  Tape  6 con- 
tains the  time  varying  gain  schedxile  used  for  the  optimal 
performance  and  used  as  input  to  compute  the  constant  gain 

The  input  data  is  listed  and  explained  in  the  program 
listing. 
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PROGRAM  GAIN (INPUT, OUTPUT, TAPE5,TA3E6) 
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